In [ ]:
import matplotlib.pyplot as plt
import cobra
from cobra.io import validate_sbml_model
import importlib
from xml.etree import ElementTree
import utils.Model_correction as mc
import sys
import utils.model_maj as mj
import utils.viz_utils as vu
#import cplex
import cbmpy
import pandas as pd
import plotly.express as px

pyparsing import

INFO: No xlwt module available, Excel spreadsheet creation disabled
CBGLPK based on swiglpk: not all methods implimented yet! 5.0

*****
Using CPLEX
*****

doFBAMinSum not available with GLPK

INFO: No xlrd module available, Excel spreadsheet reading disabled


***********************************************************************
* Welcome to CBMPy (0.8.4) - PySCeS Constraint Based Modelling        *
*                http://cbmpy.sourceforge.net                         *
* Copyright(C) Brett G. Olivier 2014 - 2020                           *
* Systems Biology Lab, Vrije Universiteit Amsterdam                   *
* Amsterdam, The Netherlands                                          *
* CBMPy is developed as part of the BeBasic MetaToolKit Project       *
* Distributed under the GNU GPL v 3.0 licence, see                    *
* LICENCE (supplied with this release) for details                    *
***********************************************************************
In [ ]:
importlib.reload(mj)
importlib.reload(vu)
importlib.reload(mc)
Out[ ]:
<module 'utils.Model_correction' from '/home/jmuller/Documents/Stage_metabo/code/utils/Model_correction.py'>

Chargement des modèles¶

In [ ]:
HepG2, errors = validate_sbml_model("../models_storage/Hep-G2_v3.xml")
In [ ]:
iHep, errors = validate_sbml_model("../models_storage/iHepatocytes2322_Cobra.xml")
In [ ]:
HepG2 = mj.maj("iHepatocytes2322_Cobra.xml", "Hep-G2_v2.xml", "../models_storage/", bounds_check = False, genes_id_copy= False, alt_gene_ids= False, metab_id_check= False, bounds_value_check= False, subsystem_copy= True )
In [ ]:
# Bidouillage pour ajouter des subsystems aux réactions d'échange. (Va ajouter "exchange reactions" normalement.)
for r in HepG2.reactions :
    if "EX_" in r.id :
        m = [m_id for m_id, _ in HepG2.reactions.get_by_id(r.id).metabolites.items()]
        sub = [r.subsystem for r in m[0].reactions if "EX_" not in r.id]
        r.subsystem = sub[0]

Bounds modifications¶

In [ ]:
HepG2.reactions.EX_m01965x.bounds = (-1000.0,-1000.0) #Glucose exchange
In [ ]:
HepG2.reactions.HMR_3883.bounds = (0.0,0.0) #Serine --> Pyruvate exchange
HepG2.reactions.EX_m02630x.bounds = (-1000.0,1000.0) #O2 exchange
HepG2.reactions.EX_m02819x.bounds = (0.0,1000.0) #Pyruvate exchange
HepG2.reactions.EX_m01910x.bounds = (0.0,0.0) #Galactose exchange
HepG2.reactions.EX_m01743x.bounds = (0.0, 1000.0) #D-Ribulose exchange
HepG2.reactions.EX_m01962x.bounds = (0.0,1000.0) #glucosamine exchange
HepG2.reactions.HMR_4316.bounds = (-1000.0,1000.0) # Glucose --> D-Glucitol
HepG2.reactions.EX_m02896x.bounds = (0.0,1000.0) # Serine intake
HepG2.reactions.EX_m01682x.bounds = (-1000.0,0.0) # Glucitol secretion blocked. Kept the intake just in case.
HepG2.reactions.EX_m01840x.bounds = (-1000.0,0.0) #Fructose exchange
#HepG2.reactions.HMR_4930.bounds = (-1000.0,1000.0) # Pyruvate transfer from cytoplasm to peroxysome
#HepG2.reactions.HMR_4281.bounds = (0.0,0.0) # Fermentation in peroxysome
In [ ]:
#HepG2.reactions.EX_m02403x.bounds = (0.0,1000.0) #Lactate exchange
HepG2.reactions.EX_m02819x.bounds = (0.0,1000.0) #Pyruvate exchange
#HepG2.reactions.EX_m01840x.bounds = (0.0,0.0) #Fructose exchange
In [ ]:
HepG2.reactions.HMR_4281.bounds = (-1000.0,0.0) # Pyruvate fermentation in peroxysome
In [ ]:
HepG2.reactions.EX_m02630x.bounds = (0.0,1000.0) #O2 exchange
In [ ]:
HepG2.reactions.EX_m02630x.bounds = (-1000.0,1000.0) #O2 exchange
HepG2.reactions.EX_m02819x.bounds = (0.0,1000.0) #Pyruvate exchange
HepG2.reactions.EX_m01910x.bounds = (0.0,0.0) #Galactose exchange
HepG2.reactions.EX_m01743x.bounds = (0.0, 1000.0) #D-Ribulose exchange
HepG2.reactions.EX_m01962x.bounds = (0.0,1000.0) #glucosamine exchange
#HepG2.reactions.EX_m01965x.bounds = (-1000,1000.0) #Glucose exchange
#HepG2.reactions.EX_m01286x.bounds = (0.0,0.0) #ADP-Glucose exchange
#HepG2.reactions.EX_m01840x.bounds = (-1000.0,0.0) #Fructose exchange
#HepG2.reactions.EX_m01743x.bounds = (-1000.0,1000.0) #D-Ribulose exchange
#HepG2.reactions.EX_m02453x.bounds = (0.0,0.0) #Mannose exchange
#HepG2.reactions.EX_m02470x.bounds = (0.0,0.0) #Methanol exchange

# --> After Methanol cut, Glucose uptake flux went from 0.0 to -190.0
# Without O2, glucose uptake went down to -23

Optimization¶

In [ ]:
HepG2.objective = "biomass_components"
sol = HepG2.optimize()
sol
Out[ ]:
Optimal solution with objective value 27.150
fluxes reduced_costs
EX_m00097x -0.000000 0.000000e+00
EX_m00157x 952.717900 0.000000e+00
EX_m00228x 0.000000 0.000000e+00
EX_m00242x 0.000000 0.000000e+00
EX_m00266x 0.000000 0.000000e+00
... ... ...
HMR_9715 0.000000 0.000000e+00
HMR_9730 0.000000 0.000000e+00
HMR_9736 36.365246 0.000000e+00
biomass_components 27.150400 4.163336e-17
artificial_biomass 0.000000 -2.000000e+00

5906 rows × 2 columns

In [ ]:
new_iHep.objective = "biomass_components"
sol = new_iHep.optimize()
sol

See exchanges¶

In [ ]:
[(m.name, m.id) for m in HepG2.reactions.biomass_components.metabolites]
In [ ]:
importlib.reload(vu)
vu.print_exchanges(HepG2, filter = "all")

run_parcours¶

In [ ]:
f = vu.run_parcours(HepG2.reactions.EX_m01965x)
print(len(f))
#m01682c glucitol

see parcours¶

In [ ]:
for reaction, flux in f.items() :
    vu.print_reactions(HepG2.reactions.get_by_id(reaction), flux)
    print(f"FLUX : {flux} --- ID : {HepG2.reactions.get_by_id(reaction).id} --- compartment : {HepG2.reactions.get_by_id(reaction).compartments} \n\n---\n\n")
In [ ]:
vu.print_reactions(HepG2.reactions.biomass_components)

FVA¶

In [ ]:
from cobra.flux_analysis import flux_variability_analysis
In [ ]:
fva = flux_variability_analysis(HepG2)
In [ ]:
indispensable = []
unisens = []
reversible = []
unknown = []


for reaction_id, flux_values in fva.iterrows() :
    min_flux = flux_values["minimum"]
    max_flux = flux_values["maximum"]
    difference = abs(max_flux - min_flux)

    if max_flux >0.0 and min_flux > 0.0 :
        f = 1.0
    elif max_flux > 0.0 and min_flux < 0.0 :
        f = 0.0
    elif max_flux < 0.0 and min_flux < 0.0 :
        f = -1.0
    else :
        f = 0.0


    if difference < 0.01 and min_flux + max_flux != 0.0 and abs(min_flux + max_flux) > 0.01 :
       indispensable.append((reaction_id, flux_values))

    elif min_flux < 0.0 and max_flux > 0.0 :
        reversible.append((reaction_id, flux_values))

    elif min_flux > 0.0 and max_flux > 0.0 or min_flux < 0.0 and max_flux < 0.0 :
        if difference > 0.01 :
            unisens.append((reaction_id, flux_values))
    else :
        unknown.append((reaction_id, flux_values))

print("\n--------INDISPENSABLE--------\n") 
for reaction_id, flux_values in indispensable : 
    min_flux = flux_values["minimum"]
    max_flux = flux_values["maximum"]   
    vu.print_reactions(HepG2.reactions.get_by_id(reaction_id), f)
    print(f"\nMIN FLUX -- {min_flux} === {max_flux} -- MAX FLUX\n\n")

print("\n--------UNISENS--------\n")
for reaction_id, flux_values in unisens :    
    min_flux = flux_values["minimum"]
    max_flux = flux_values["maximum"]
    vu.print_reactions(HepG2.reactions.get_by_id(reaction_id), f)
    print(f"\nMIN FLUX -- {min_flux} === {max_flux} -- MAX FLUX\n\n")

print("\n--------REVERSIBLE--------\n")
for reaction_id, flux_values in reversible :  
    min_flux = flux_values["minimum"]
    max_flux = flux_values["maximum"]  
    vu.print_reactions(HepG2.reactions.get_by_id(reaction_id), f)
    print(f"\nMIN FLUX -- {min_flux} === {max_flux} -- MAX FLUX\n\n")

if len(unknown) >0 :
    print("\n--------UNKNOWN--------\n")
    for reaction_id, flux_values in unknown :    
        min_flux = flux_values["minimum"]
        max_flux = flux_values["maximum"]
        vu.print_reactions(HepG2.reactions.get_by_id(reaction_id), f)
        print(f"\nMIN FLUX -- {min_flux} === {max_flux} -- MAX FLUX\n\n")

Visualisation¶

In [ ]:
df_C, df_R, df_S, df_M, df_P, df_X, df_L, df_G, df_N = vu.build_reaction_df(HepG2)
In [ ]:
flux_filter = 0.0
In [ ]:
fig = px.treemap(df_C.loc[(df_C["flux"] >= flux_filter)] , path=['subSystem', 'id'], 
                 values='flux', color='subSystem',color_discrete_sequence=px.colors.qualitative.Light24_r)
fig.update_layout(title_text="Cytoplasm reactions", font_size=12)
fig.show(renderer='notebook')
In [ ]:
fig = px.treemap(df_R.loc[df_R["flux"] >= flux_filter], path=['subSystem', 'id'], 
                 values='flux', color='subSystem',color_discrete_sequence=px.colors.qualitative.Light24_r)
fig.update_layout(title_text="Endoplasmic Reticulum reactions", font_size=12)
fig.show(renderer='notebook')
In [ ]:
fig = px.treemap(df_S.loc[df_S["flux"] >= flux_filter], path= ['subSystem', 'id'], 
                 values='flux', color='subSystem',color_discrete_sequence=px.colors.qualitative.Light24_r)
fig.update_layout(title_text="Endomembraneous reactions", font_size=12)
fig.show(renderer='notebook')
In [ ]:
fig = px.treemap(df_M.loc[df_M["flux"] >= flux_filter], path= [px.Constant("Mitochondria"),'subSystem', 'id'], 
                 values='flux', color='subSystem',color_discrete_sequence=px.colors.qualitative.Light24_r)
fig.update_layout(title_text="Mitochondrial reactions", font_size=12)
fig.show(renderer='notebook')
In [ ]:
fig = px.treemap(df_P.loc[df_P["flux"] >= flux_filter], path= [px.Constant("Peroxysome"),'subSystem', 'id'], 
                 values='flux', color='subSystem',color_discrete_sequence=px.colors.qualitative.Light24_r)
fig.update_layout(title_text="Peroxysomal reactions", font_size=12)
fig.show(renderer='notebook')
In [ ]:
fig = px.treemap(df_X.loc[df_X["flux"] >= flux_filter], path= ['subSystem', 'id'], 
                 values='flux', color='subSystem',color_discrete_sequence=px.colors.qualitative.Light24_r)
fig.update_layout(title_text="Exchange reactions", font_size=12)
fig.show(renderer='notebook')
In [ ]: